MODULE indexx

  IMPLICIT NONE
  PRIVATE
  INTEGER, PARAMETER :: NSWITCH=60

  ! Callable routines in this module 
  PUBLIC :: indexxx

  ! Sorting routine provided for:
  ! REAL, INTEGER, and INTEGER*8
  INTERFACE indexxx
     MODULE PROCEDURE indexxxr
     MODULE PROCEDURE indexxxi
     MODULE PROCEDURE indexxxi8
  END INTERFACE

CONTAINS

! **** Sorting routines for REAL arrays ****

  SUBROUTINE  indexxxr(n,arr,indx) 
    !     
    IMPLICIT NONE
    INTEGER n,NSWITCH
    REAL arr(n)
    INTEGER indx(n)
    !     
    IF (n.GT.NSWITCH) THEN     !if n is large use the Numerical recipes routine
       CALL indexxr(n,arr,indx) !heapsort from numerical recipes 
    ELSE                       !simple shell sort which is faster for
                               !small n but scales as n^3/2 in worst case
       CALL indexshr(n,arr,indx)!rather than the n log(n) of heapsort
    END IF
    !     
    RETURN
  END SUBROUTINE indexxxr
  !     
  !     Indexing or array arr(n) using shell sort
  !     
  SUBROUTINE  indexshr(n,arr,indx)
    !     
    IMPLICIT NONE
    INTEGER n,m
    REAL arr(n),aln2i,tiny
    PARAMETER(TINY=1.0e-05,ALN2I=1.0/0.69314718)
    INTEGER indx(n),nn,lognb2,l,i,tindx,j,k
    !     
    DO i=1,n
       indx(i)=i
    END DO
    IF(n.GE.2) THEN
       lognb2=INT(LOG(float(n))*ALN2I+TINY)
       m=n
       DO nn=1,lognb2
          m=m/2
          k=n-m
          DO j=1,k
             i=j
3            CONTINUE
             l=i+m
             IF(arr(indx(l)).LT.arr(indx(i))) THEN
                tindx=indx(i)
                indx(i)=indx(l)
                indx(l)=tindx
                i=i-m
                IF (i.GE.1) GOTO 3
             END IF
          END DO
       END DO
    END IF
    !     
    RETURN
  END SUBROUTINE indexshr
  !
  !     Numerical Recipes heapsort
  SUBROUTINE INDEXXR(N,ARRIN,INDX)
    IMPLICIT NONE
    INTEGER N
    INTEGER INDX(N),J,L,I,IR,INDXT
    REAL ARRIN(N),Q
    DO J=1,N
       INDX(J)=J
    ENDDO
    IF(n.LE.1) RETURN
    L=N/2+1
    IR=N
10  CONTINUE
    IF(L.GT.1)THEN
       L=L-1
       INDXT=INDX(L)
       Q=ARRIN(INDXT)
    ELSE
       INDXT=INDX(IR)
       Q=ARRIN(INDXT)
       INDX(IR)=INDX(1)
       IR=IR-1
       IF(IR.EQ.1)THEN
          INDX(1)=INDXT
          RETURN
       ENDIF
    ENDIF
    I=L
    J=L+L
20  IF(J.LE.IR)THEN
       IF(J.LT.IR)THEN
          IF(ARRIN(INDX(J)).LT.ARRIN(INDX(J+1)))J=J+1
       ENDIF
       IF(Q.LT.ARRIN(INDX(J)))THEN
          INDX(I)=INDX(J)
          I=J
          J=J+J
       ELSE
          J=IR+1
       ENDIF
       GO TO 20
    ENDIF
    INDX(I)=INDXT
    GO TO 10
  END SUBROUTINE INDEXXR

!
! **** Sorting routines for INTEGER arrays ****
!     
  SUBROUTINE  indexxxi(n,arr,indx) 
    !     
    IMPLICIT NONE
    INTEGER n
    INTEGER arr(n)
    INTEGER indx(n)
    !     
    IF (n.GT.NSWITCH) THEN      !if n is large use the Numerical recipes
       !routine
       CALL indexxi(n,arr,indx) !heapsort from numerical recipes 
    ELSE                        !simple shell sort which is faster
       !for small n but scales as n^3/2 in worst case
       CALL indexshi(n,arr,indx) !rather than the n log(n) of heapsort
    END IF
    !     
    RETURN
  END SUBROUTINE indexxxi
  !     
  !     Indexing or array arr(n) using shell sort
  !     
  SUBROUTINE  indexshi(n,arr,indx)
    !     
    IMPLICIT NONE
    INTEGER n,m
    INTEGER arr(n)
    REAL aln2i,tiny
    PARAMETER(TINY=1.0e-05,ALN2I=1.0/0.69314718)
    INTEGER indx(n),nn,lognb2,l,i,tindx,j,k
    !     
    DO i=1,n
       indx(i)=i
    END DO
    IF(n.GE.2) THEN
       lognb2=INT(LOG(float(n))*ALN2I+TINY)
       m=n
       DO nn=1,lognb2
          m=m/2
          k=n-m
          DO j=1,k
             i=j
3            CONTINUE
             l=i+m
             IF(arr(indx(l)).LT.arr(indx(i))) THEN
                tindx=indx(i)
                indx(i)=indx(l)
                indx(l)=tindx
                i=i-m
                IF (i.GE.1) GOTO 3
             END IF
          END DO
       END DO
    END IF
    !     
    RETURN
  END SUBROUTINE indexshi
  !     
  !     Numerical Recipes heapsort
  SUBROUTINE INDEXXI(N,ARRIN,INDX)
    IMPLICIT NONE
    INTEGER N
    INTEGER INDX(N),J,L,I,IR,INDXT
    INTEGER ARRIN(N),Q
    DO J=1,N
       INDX(J)=J
    ENDDO
    IF(n.LE.1) RETURN
    L=N/2+1
    IR=N
10  CONTINUE
    IF(L.GT.1)THEN
       L=L-1
       INDXT=INDX(L)
       Q=ARRIN(INDXT)
    ELSE
       INDXT=INDX(IR)
       Q=ARRIN(INDXT)
       INDX(IR)=INDX(1)
       IR=IR-1
       IF(IR.EQ.1)THEN
          INDX(1)=INDXT
          RETURN
       ENDIF
    ENDIF
    I=L
    J=L+L
20  IF(J.LE.IR)THEN
       IF(J.LT.IR)THEN
          IF(ARRIN(INDX(J)).LT.ARRIN(INDX(J+1)))J=J+1
       ENDIF
       IF(Q.LT.ARRIN(INDX(J)))THEN
          INDX(I)=INDX(J)
          I=J
          J=J+J
       ELSE
          J=IR+1
       ENDIF
       GO TO 20
    ENDIF
    INDX(I)=INDXT
    GO TO 10
  END SUBROUTINE INDEXXI

! **** Sorting routines for INTEGER*8 arrays ****

  SUBROUTINE  indexxxi8(n,arr,indx) 
    !     
    IMPLICIT NONE
    INTEGER n
    INTEGER*8 arr(n)
    INTEGER indx(n)
    !     
    IF (n.GT.NSWITCH) THEN    
       CALL indexxi8(n,arr,indx) 
    ELSE                      
       CALL indexshi8(n,arr,indx) 
    END IF
    !     
    RETURN
  END SUBROUTINE indexxxi8
  !     
  !     Indexing or array arr(n) using shell sort
  !     
  SUBROUTINE  indexshi8(n,arr,indx)
    !     
    IMPLICIT NONE
    INTEGER n,m
    INTEGER*8 arr(n)
    REAL aln2i,tiny
    PARAMETER(TINY=1.0e-05,ALN2I=1.0/0.69314718)
    INTEGER indx(n),nn,lognb2,l,i,tindx,j,k
    !     
    DO i=1,n
       indx(i)=i
    END DO
    IF(n.GE.2) THEN
       lognb2=INT(LOG(float(n))*ALN2I+TINY)
       m=n
       DO nn=1,lognb2
          m=m/2
          k=n-m
          DO j=1,k
             i=j
3            CONTINUE
             l=i+m
             IF(arr(indx(l)).LT.arr(indx(i))) THEN
                tindx=indx(i)
                indx(i)=indx(l)
                indx(l)=tindx
                i=i-m
                IF (i.GE.1) GOTO 3
             END IF
          END DO
       END DO
    END IF
    !     
    RETURN
  END SUBROUTINE indexshi8
  !     

  !
  SUBROUTINE INDEXXI8(N,ARRIN,INDX)
    IMPLICIT NONE
    INTEGER N
    INTEGER INDX(N),J,L,I,IR,INDXT
    INTEGER*8 ARRIN(N)
    INTEGER*8 Q
    DO J=1,N
       INDX(J)=J
    ENDDO
    IF(n.LE.1) RETURN
    L=N/2+1
    IR=N
10  CONTINUE
    IF(L.GT.1)THEN
       L=L-1
       INDXT=INDX(L)
       Q=ARRIN(INDXT)
    ELSE
       INDXT=INDX(IR)
       Q=ARRIN(INDXT)
       INDX(IR)=INDX(1)
       IR=IR-1
       IF(IR.EQ.1)THEN
          INDX(1)=INDXT
          RETURN
       ENDIF
    ENDIF
    I=L
    J=L+L
20  IF(J.LE.IR)THEN
       IF(J.LT.IR)THEN
          IF(ARRIN(INDX(J)).LT.ARRIN(INDX(J+1)))J=J+1
       ENDIF
       IF(Q.LT.ARRIN(INDX(J)))THEN
          INDX(I)=INDX(J)
          I=J
          J=J+J
       ELSE
          J=IR+1
       ENDIF
       GO TO 20
    ENDIF
    INDX(I)=INDXT
    GO TO 10
  END SUBROUTINE INDEXXI8

END MODULE indexx
